% 本例程用于在第二节说明不同采样率下包络起伏情况对比

clear;
close all;

load('non_smoth_weights.mat', 'non_smoth_weights');
% load('smoth_weights.mat', 'smoth_weights');
load('smoth_weights_pca_l1.mat', 'smoth_weights');

weightLen = length(non_smoth_weights);
os2 = 2;
os4 = 4;

ampNonsmothOs1 = abs(ifft(non_smoth_weights));
ampNonsmothOs2 = 2*abs(ifft(non_smoth_weights, os2*weightLen));
ampNonsmothOs4 = 4*abs(ifft(non_smoth_weights, os4*weightLen));

angNonsmothOs1 = angle(ifft(non_smoth_weights));
angNonsmothOs2 = angle(ifft(non_smoth_weights, os2*weightLen));
angNonsmothOs4 = angle(ifft(non_smoth_weights, os4*weightLen));

unwrapedAngNonsmothOs1 = unwrap_increasing(angNonsmothOs1);
unwrapedAngNonsmothOs2 = unwrap_increasing(angNonsmothOs2);
unwrapedAngNonsmothOs4 = unwrap_increasing(angNonsmothOs4);

angDiffNonsmothOs1 = unwrapedAngNonsmothOs1(2:end) - unwrapedAngNonsmothOs1(1:end-1);
angDiffNonsmothOs2 = unwrapedAngNonsmothOs2(2:end) - unwrapedAngNonsmothOs2(1:end-1);
angDiffNonsmothOs4 = unwrapedAngNonsmothOs4(2:end) - unwrapedAngNonsmothOs4(1:end-1);


ampSmothOs1 = abs(ifft(smoth_weights));
ampSmothOs2 = 2*abs(ifft(smoth_weights, os2*weightLen));
ampSmothOs4 = 4*abs(ifft(smoth_weights, os4*weightLen));

angSmothOs1 = angle(ifft(smoth_weights));
angSmothOs2 = angle(ifft(smoth_weights, os2*weightLen));
angSmothOs4 = angle(ifft(smoth_weights, os4*weightLen));

unwrapedAngSmothOs1 = unwrap_increasing(angSmothOs1);
unwrapedAngSmothOs2 = unwrap_increasing(angSmothOs2);
unwrapedAngSmothOs4 = unwrap_increasing(angSmothOs4);

angDiffSmothOs1 = unwrapedAngSmothOs1(2:end) - unwrapedAngSmothOs1(1:end-1);
angDiffSmothOs2 = unwrapedAngSmothOs2(2:end) - unwrapedAngSmothOs2(1:end-1);
angDiffSmothOs4 = unwrapedAngSmothOs4(2:end) - unwrapedAngSmothOs4(1:end-1);

index1 = (1:os4:os4*weightLen);
index2 = (1:os2:os4*weightLen);
index4 = (1:1:os4*weightLen);

figure(1);hold on;
title("包络对比")
subplot(231)
h = stem(index1, ampNonsmothOs1);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');

subplot(234)
h = stem(index1, ampSmothOs1);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');

subplot(232)
h = stem(index2, ampNonsmothOs2);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');

subplot(235)
h = stem(index2, ampSmothOs2);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');

subplot(233)
h = stem(index4, ampNonsmothOs4);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');

subplot(236)
h = stem(index4, ampSmothOs4);
h.Marker = 'o';                  % 点的形状为圆
h.MarkerFaceColor = 'b';         % 实心颜色
h.MarkerEdgeColor = 'b';         % 边框颜色
h.MarkerSize = 2.5;                % 控制点的大小（直径，单位是像素）
ylim([0 2])
xlim([49 112])
ytickformat('%.1f');


figure(2);hold on;
title("相邻元素的相位变化对比")
subplot(231)
h = plot(index1, angNonsmothOs1);
ylim([-pi pi])

subplot(234)
h = plot(index1, angSmothOs1);
ylim([-pi pi])
ytickformat('%.1f');

subplot(232)
h = plot(index2, angNonsmothOs2);
ylim([-pi pi])

subplot(235)
h = plot(index2, angSmothOs2);
ylim([-pi pi])
ytickformat('%.1f');

subplot(233)
h = plot(index4, angNonsmothOs4);
ylim([-pi pi])
ytickformat('%.1f');

subplot(236)
h = plot(index4, angSmothOs4);
ylim([-pi pi])
ytickformat('%.1f');


figure(3);hold on;
title("相邻元素的相位一阶差分对比")
subplot(231);hold on;
h = plot(index1(1:end-1), angDiffNonsmothOs1);
h = plot(index2(1:end-1), 2*angDiffNonsmothOs2);
h = plot(index4(1:end-1), 4*angDiffNonsmothOs4);
% ylim([0 2*pi])

subplot(234);hold on;
h = plot(index1(1:end-1), angDiffSmothOs1);
h = plot(index2(1:end-1), 2*angDiffSmothOs2);
h = plot(index4(1:end-1), 4*angDiffSmothOs4);
ylim([0 2*pi])
ytickformat('%.1f');

subplot(232)
h = plot(index2(1:end-1), angDiffNonsmothOs2);
ylim([0 2*pi])

subplot(235)
h = plot(index2(1:end-1), angDiffSmothOs2);
ylim([0 2*pi])
ytickformat('%.1f');

subplot(233)
h = plot(index4(1:end-1), angDiffNonsmothOs4);
ylim([0 2*pi])
ytickformat('%.1f');

subplot(236)
h = plot(index4(1:end-1), angDiffSmothOs4);
ylim([0 2*pi])
ytickformat('%.1f');


% % 自相关研究
nonSmothOs1 = ifft(non_smoth_weights);
nonSmothOs2 = ifft(non_smoth_weights, os2 * weightLen);
nonSmothOs4 = ifft(non_smoth_weights, os4 * weightLen);

acfNonsmothOs1 = conv(conj(flip(nonSmothOs1)), nonSmothOs1);
acfNonsmothOs2 = conv(conj(flip(nonSmothOs2)), nonSmothOs2);
acfNonsmothOs4 = conv(conj(flip(nonSmothOs4)), nonSmothOs4);

acfIndex1 = 4-4*weightLen:4:4*weightLen-4;
acfIndex2 = 2-4*weightLen:2:4*weightLen-2;
acfIndex4 = 1-4*weightLen:1:4*weightLen-1;

figure(4);hold on;legend();
title("非平滑序列不同采样率下ACF");
plot(acfIndex1, 20*log10(abs(acfNonsmothOs1) / max(abs(acfNonsmothOs1))), 'DisplayName', 'Os=1时的ACF');
plot(acfIndex2, 20*log10(abs(acfNonsmothOs2) / max(abs(acfNonsmothOs2))), 'DisplayName', 'Os=2时的ACF');
plot(acfIndex4, 20*log10(abs(acfNonsmothOs4) / max(abs(acfNonsmothOs4))), 'DisplayName', 'Os=4时的ACF');


% % 平滑恒模序列自相关研究
SmothOs1 = ifft(smoth_weights);
SmothOs2 = ifft(smoth_weights, os2 * weightLen);
SmothOs4 = ifft(smoth_weights, os4 * weightLen);

acfSmothOs1 = conv(conj(flip(SmothOs1)), SmothOs1);
acfSmothOs2 = conv(conj(flip(SmothOs2)), SmothOs2);
acfSmothOs4 = conv(conj(flip(SmothOs4)), SmothOs4);

acfIndex1 = 4-4*weightLen:4:4*weightLen-4;
acfIndex2 = 2-4*weightLen:2:4*weightLen-2;
acfIndex4 = 1-4*weightLen:1:4*weightLen-1;

figure(5);hold on;legend();
title("平滑序列不同采样率下ACF");
plot(acfIndex1, 20*log10(abs(acfSmothOs1) / max(abs(acfSmothOs1))), 'DisplayName', 'Os=1时的ACF');
plot(acfIndex2, 20*log10(abs(acfSmothOs2) / max(abs(acfSmothOs2))), 'DisplayName', 'Os=2时的ACF');
plot(acfIndex4, 20*log10(abs(acfSmothOs4) / max(abs(acfSmothOs4))), 'DisplayName', 'Os=4时的ACF');






function y = unwrap_increasing(phase)
    y = zeros(size(phase));
    y(1) = phase(1);
    for k = 2:length(phase)
        delta = phase(k) - phase(k-1);
        while delta < 0
            phase(k) = phase(k) + 2*pi;
            delta = phase(k) - phase(k-1);
        end
        y(k) = phase(k);
    end
end

